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ABSTRACT 


We  have  designed  a  cubic  spline  wavelet  decomposition  for  the  Sobolev  space  Hq(I)  where 
/  is  a  bounded  interval.  Based  on  a  special  “point-wise  orthogonality”  of  the  wavelet  basis 
functions,  a  fast  Discrete  Wavelet  Transform  (DWT)  is  constructed.  This  DV  f  transform 
will  map  discrete  samples  of  a  function  to  its  wavelet  expansion  coefficients  in  O(A'logiV) 
operations.  Using  this  transform,  we  propose  a  collocation  method  for  the  initial  value 
boundary  problem  of  nonlinear  PDE’s.  Then,  we  test  the  efficiency  of  the  DWT  transform 
and  apply  the  collocation  method  to  solve  linear  and  nonlinear  PDE’s. 
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1  Introduction 


Wavelet  approximations  have  attracted  much  attention  as  a  potential  efficient  numerical 
technique  for  the  solutions  of  partial  differential  equations  [1]  -  [6].  Because  of  their  advan¬ 
tageous  properties  of  localizations  in  both  space  and  frequency  domains  [8]  -  [10],  wavelets 
seem  to  be  a  great  candidate  for  adaptive  schemes  for  solutions  which  vary  dramatically 
both  in  space  and  time  and  develop  singularities.  However,  in  order  to  take  advantage  cf  the 
nice  properties  of  wavelet  approximations,  we  have  to  find  an  efficient  way  to  deal  with  the 
nonlinearity  and  general  boundary  conditions  in  the  PDE’s.  After  all,  most  of  the  problems 
of  fluid  dynamics  and  electromagnetics,  which  involve  solutions  with  quite  different  scales, 
are  governed  by  nonlinear  PDE’s  with  complicated  boundary  conditions.  Therefore,  it  is  our 
objective  here  to  address  these  issues  when  designing  wavelet  approximations  and  numerical 
schemes  for  nonlinear  PDE’s. 

We  will  present  a  new  wavelet  collocation  method  designed  to  solve  nonlinear  time  evo¬ 
lution  problems.  The  key  component  in  this  collocation  method  is  a  so-called  “Discrete 
Wavelet  Transform”  (DWT)  which  maps  a  solution  between  the  physical  space  and  the 
wavelet  coefficient  space.  The  wavelet  decomposition  is  based  on  a  new  cubic  spline  wavelet 
for  Hq{I)  where  /  is  a  bounded  interval  [11].  In  order  to  treat  the  boundary  conditions  an  ex¬ 
tra  boundary  scaling  function  4>b{x)  and  a  boundary  wavelet  i/>t(x)  have  been  used.  A  special 
“pointwise  orthogonality”  (  see(3.7))  of  the  wavelet  functions  ipj,k(z)  results  in  O(NlogN) 
operations  for  the  DWT  transform  where  N  is  the  total  number  of  unkncwns.  Therefore,  the 
nonlinear  term  in  the  PDE  can  be  easily  treated  in  the  physical  space,  and  the  derivatives  of 
those  nonlinear  terms  then  computed  in  the  wavelet  space.  As  a  result,  collocation  methods 
will  provide  the  flexibility  of  handling  nonlinearity  (and,  also  the  implementation  of  various 
boundary  conditions)  which  usually  are  not  shared  by  Galerkin  type  wavelet  methods  and 
finite  element  methods. 

The  rest  of  this  paper  is  divided  into  the  following  five  sections.  In  section  2,  we  introduce 
the  cubic  scaling  functions  ^(x),  <f>b(x)  and  their  wavelet  functions  */>(x),  */’<,(. r).  A  multires- 


olution  analysis  (MR A)  and  its  corresponding  wavelet  decomposition  of  the  Sobolev  space 
//o  ( / )  are  constructed  using  <f>b(x)  and  t/’(x),  V’i>(j"  )•  Then,  we  show  how-  to  construct 

a  wavelet  approximation  for  function  in  Sobolev  space  // 2 ( / )  which  the  solutions  of  PDF's 
will  belong  to.  In  Section  3,  we  discuss  the  fast  discrete  wavelet  transform  (DWT)  between 
functions  and  their  wavelet  coefficients.  In  Section  4,  we  discuss  the  derivative  matrix  V  asso¬ 
ciated  with  wavelet  interpolations.  In  Section  5,  we  present  the  wavelet  collocation  methods 
for  nonlinear  time  evolution  PDE’s.  In  Section  6,  we  give  the  CPU  time  performance  of  the 
DWT  transforms  and  the  numerical  results  of  the  wavelet  collocation  methods  for  linear  and 
nonlinear  PDE’s,  and  a  conclusion  is  given  in  Section  7. 

2  Scaling  functions  </>(z),  <j>b{x)  and  wavelet  functions 

Let  /  denote  any  finite  interval,  say  /  =  [0,  L]  and  L  is  a  positive  integer  (for  the  sake  of 
simplicity,  we  assume  that  L  >  4),  and  H2(I)  and  Hq(I)  denote  the  following  two  Sobolev 
spaces  with  finite  L 2  norm  for  up  to  the  second  derivatives,  i.e. 

H‘(l)  =  {/(*),*  €  1\  ||/|i)||,  <  oo,i  =  0, 1,2}  (2.1) 

Hun  =  {/(*>  e  h\i)  i  m  =  m  =  m  =  rw  =  •>}.  (2.2) 

It  can  be  easily  checked  [7]  that  Hq(I)  is  a  Hilbert  space  with  the  inner  product 

<  f,9  >=  f"{x)g"(x)  dx,  (2.3) 

thus, 

lll/lll  =  v/</./>  (2.4) 

provides  a  norm  for  Hq(I). 

In  order  to  generate  a  multiresolution  for  Sobolev  space  Hq(I),  we  consider  two  scaling 
functions,  an  interior  scaling  function  <f>(x)  and  a  boundary  scaling  function  <f>b{x)  (see  Figure 
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I) 


«U)  =  .iV<(i)  =  i^:(j)(-lV(x-;)3+ 

0  J=0  V  ' 

Mx)  =  ^4  -  ^  -  1)+  -  j(-r  -  2)+ 

where  A^j*)  is  the  4th  order  B-spline  [13]  and  for  any  real  number  n 

„  _  |  r"  if  x  >  0 
+  1  0  otherwise. 

In  a  pair  they  satisfy  the  following  two-scale  relationship, 


Lemma  1 


o(x)  =  X]23(j)d>( 2x-k) 


=  3--[<Pb(2x)  +  ^2  ~  k) 

h~0 

,  ,,  3  ..  3  17  13 

here  4_i  =  -,/Jo  =  =  —  ./h  =  — 

4  b  4  4 


We  summarize  some  properties  of  0(x)  and  Ob(x)  in  the  following  lemma. 
Lemma  2  Let  0(x)  and  Ob{x)  be  defined  as  in  (2.5)  and  (2.6),  then  we  have 


(1)  supp{(p(x))  =  [0,4]; 

(2)  supp(<t>h(x))  =  [0,3]; 

(3)  <j>(x),M*)€  //02(/); 

(4)  0'(1)  =  -0'( 3)  =  i  0'(2)  =  0,  0'6(1)  =  ^,0'6(2)  =  -i; 

(5)  0(1)  =  0(3)  =  i  0(2)  =  06(1)  =  ^  0fc(2)  =  I 

6  3  12  6 


(2.10) 

(2.11) 

(2.12) 


For  any  j,  fc  G  Z,  we  define 


0j,*(x)  =  0(2Jx  -  A:).  06,7(-r)  =  06(2Jx). 


(2.13) 


And  for  each  j,  let  V3  be  the  closure  under  norm  j 
{ dj'ic(x),  0  <  k  <  2J  L  —  4,  <t>b,j(x).  6b,j[L  —  x)},  namely 


in  (2.4)  of  the  linear  span  of 


V}  =  span{3Jjc(x),  0  <  k  <  2 JL  -  4,  d>b,}(x),  <pbj{L  ~  x)}. 


(2.14) 
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Theorem  1  Let  V ).]  £  Z+  he  the  linear  span  of  (2.14),  then  V)  forms  a  mult iresolution 
analysis  (MRA)  for  H^(I)  equipped  with  norm  (2.4)  in  the  following  sense 

(i)  V0  CL.C  V2  C  ■  •  ■ ; 

(ii)  do»H2(\JjeZ+  Vj)  =  H*(I); 

(iii)  flje z+  Vj  =  lo;  and 

(iv)  for  each  j,  { ^^(x)  =  <t>('2J .v  —  k),6b,j(x)  =  Ob{‘iJ  x).<fb.j(L  —  J") }  is  an  unconditional 
basis  of  V,. 

Proof.  The  proof  for  (iii)  ami  (iv)  is  straightforward  and  omitted  here.  The  proof  for 
(i)  follows  from  (2.7)  in  Lemma  1.  In  order  to  prove  (ii),  we  recall  a  familiar  result  on 
interpolation  cubic  spline  approximation  for  smooth  functions  taken  from  [14]  and  rewritten 
for  the  proof  of  our  theorem. 

□ 

Lemma  3  Let  ir  be  the  partition  given  by  r,  =  ih,  0  <  i  <  n,h  =  and  s(.r)  be  the  cubic 
spline  interpolating  f(x)  £  C4[a,6]  at  all  points  in  7 r, 

s{ *,)  =  f{xi),  U  <  i  <  n, 

and  satisfying  the  following  boundary  conditions: 

•-'(«)  =  /'(«),  »'(b)  =  f(b).  (2.15) 

Then  s(x)  uniquely  exists  and 

ll*(r)  -  /(r) ||L*  <  r  =  0,  1,2,3  (2.1b) 

where  e0  =  =  ^,e-2  =  |,(:3  =  l. 

Proof  of  (ii)  of  Theorem  /  Let  h  =  jj,  a  —  0  and  b  =  L.  Consider  /(./•)  £  Co'(0.  L),  Since 
<7o~(0,  L)  C  C4[0,  L)  fl  H o(0,  L ),  by  Lemma  3,  there  is  an  unique  cubic  spline  corresponding 
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the  partition  ir  interpolating  /(x).  From  the  fact  that  /( 0)  =  f{L)  =  f'{ 0)  =  f'{L)  =  0,  we 
have  s(x)  in  V:  and  then 


L‘-4 

s(x)  =  C_id>6iJ(x)  +  Yi  Ck<t>j,k{x)  +  CL-3<t>b,]{L  -  x) 

(2.17) 

k= 0 

such  that 

s(xi)  =  /(x,),  0  <  i  <2 } L 

(2.18) 

where  V  =  23  L,  x,  = 

Finally,  from  (2.16)  in  Lemma  3  with  r  =  2  we  have 

l!M-/lll  =  IM|2|-/w!l<<2ll/|4,ll/22j- 

Therefore,  as  j  — >  oc,  |||s  -  /j||  — ♦  0.  This  proves  that  Q°(0,  L)  C  closH j(U;eZ+  V;)- 
Then,  Theorem  1  (ii)  follows  from  the  fact  that  (7^(0,  L)  is  dense  in  Hq(0,L). 

□ 

To  construct  a  wavelet  decomposition  of  Sobolev  space  Hq(I)  under  the  inner  product 
(2.3),  we  consider  the  following  two  wavelet  functions  (see  Figure  2)  , 

0(x)  =  2x)  +  y<£(2x  -  1)  -  ^<t>{2x  -  2)  €  V\  (2-19) 

M*)  =  ^6(2x)  - -|d>(2x)  G  v,.  (2.20) 

It  can  be  verified  that 

i/>(n)  =  =  0,  for  all  n  €  Z.  (2.21) 

Equation  (2.21)  will  be  important  in  the  construction  of  the  fast  DWT  transform  later.  And, 
equations  (2.19)  and  (2.20)  imply  that  r!>{x)  and  ifo{x)  both  belong  to  Vj.  As  usual,  we  define 
the  dilation  and  translation  of  these  two  functions 

&j,k(x)  ~  H'Fx  -  k),  j  >  0,  k  =  0,  •  ■  • ,  n}  -  3,  (2.22) 
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I.'*  ,(x)  =  t>fc( 2Jr),  (,v,f X)  =  I.\{2J{L  -  ./•))  (2.23) 

where  n,  =  1J  L.  For  the  sake  of  simplicity,  we  will  adopt  the  following;  notations 

'/V.-i(-r)  =  '-ijlx).  V>,-2  =  I-'').  (2.24 ) 

So  when  k  —  —  l,»j  —  2,  will  denote  the  two  boundary  wavelet  functions,  not  t lit* 

usual  translation  and  dilation  of  t/’(.r). 

Finally,  for  each  j  >  0,  we  define 

Wj  =  closH2  <  i /Vit(jr),  A-  =  - 1 ,  •  •  • ,  »,  -  2  >  .  (2.25) 


Theorem  2  The  IV; ,  j  >  0  defined  in  (2.25)  is  the  orthogonal  compliment  of  V)  in  Vj+1 
under  the  inner  product  (2.3),  i.e. 

(1)  VJ+x  =  Vj  ©  Wj  for  j  G  Z+.  Here  0  stands  for  14 J_ W3  under  the  inner  product  (2.3) 
and  Vj+i  =  V3  +  Wj.  Therefore, 

(2)  WjlWj+jJ  €  Z+; 

(3)  Hl(!)  =  V/o©j€z+  Wj. 

Proof.  ( 1 )  We  only  have  to  prove  Vj  ®  Wj  for  j  =  0,  namely,  for  0  <  /  <  L  —  4,  0  <  k  <  L  —  3, 


<  4>{x  —  /),  »/>(x  —  A;)  >  =  0 

(2.26) 

<  <j>{x  -  l),ipb(x)  >  =  0 

(2.27) 

<  f>b(x),xk(x  —  k)  >  =  0 

(2.28) 

<  <t>b(x),  ifb(x)  >  =  0. 

(2.29) 

Integrating  by  parts  twice  in  (2.26)  and  using  the  fact  that  i l'(x),<j>(x)  G  //(((/)  ,  we  have 

<  <f>(x  —  l),rj>(x  —  A;)  >  =  f  <p"{x  —  /)(/>"( x  —  k)  dx 

Jo 

=  <f>"{x  -  /)»//(*  -  A’)|q  -  fL  <t>V\x  -  l)r\x  -  k)  dx 

Jo 

=  —  /  d>^{x  —  l)ib'(x  —  k)  dx 

Jo 

=  -f(3)(x  -  l)i/’(x  -  A’)|q  +  [ 1  (j)i4)(x  -  Mix  —  A1 )  dx 

Jo 

—  f  <p^(. r  —  l)4'(x  —  k)  dx. 

Jo 
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Using  equation  (2.21)  and  the  identity 


o(4)ui  =  (j)  (-1  ysix-j) 

4  J=0  '  ' 

where  S(x)  is  the  Dirac-£  function,  so  we  have 

<  o{x  -  /).  v(x  -  k)  >=  ]  0 )  ( —  1  -  (k  ~  0)  =  °- 

4  j=o  v  ' 

Equations  (2.27)  -  (2.29)  can  be  shown  similarly  to  be  true  .  So  ( 1 )  follows  from  (2.19) 
and  (2.20)  and  the  fact  that  dim  l)  =  2 J  L  —3  and  dim  llj  -  2 J  L  and  dimV]+\  =  2J+1  L  —  3  = 
(2J  L  —  3)  +  2J  L  =  dim  Vj  +  dimWj ; 

(2)  follows  from  (1); 

(3)  follows  directly  from  Theorem  1  (ii). 

□ 

As  a  consequence  of  Theorem  2.  any  function  /(.r)  £  H*{I)  can  be  approximated  as 
closely  as  possible  by  a  function  fj(.r)  €  Vj  =  l'o  £•  H’o  G*  W i  W j  for  a  sufficiently  large 
j,  and  fj(x)  has  an  unique  orthogonal  decomposition 

fjix)  —  fo  +  <}o  +  <7i  +  ’  ’  ’  +  9]  (2.30) 

where  f0  G  V0,gi  G  Wj,0  <  i  <  j. 
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Approximation  for  function  in 

11  V) 

Consider  the  following  two  fun< 

t  ions 

viU)  = 

~x+  - 

+  ( 2d  i 

'/>(■/■)  = 

(1 

( 2,d2 

For  any  function  /(.r)  t  fl2(I). 

by  t  he 

Soho 

lev  embedding  t  lieorem  we  have  /  ( ,r )  t  C 1  (  / 

and.  therefore,  we  can  deline  the  following  houndary  interpolation  l,  ,/i.r).  j  >  I) 


h.jf(-r)  -  <>\>I\V-  ■>')  +  n i'lA'-'-i’)  +  «>-,//,  (2y( /-  -  •/•))  +  o,//j(2'|7.  -  .cl)  (2.ddi 

such  that 

W(0)  =  /(U).  IkJ(A)  =  /(A)  (2.dF) 

h,f'(0)  =  /'(O).  I  h,f'(L)  =  f'(L).  (2.:r>) 

It  can  he  easily  verified  that,  in  order  to  have  lh  j  f  satisfy  conditions  (2.11)  -  (2. do),  we 
have  to  take 

«1  =  -  ^/(0)-  ^  =  /(0)  (2.d(i) 

fU)  d 

°3  =  -T^r~-/(^).  n4  =  f(L). 


In  many  situations  we  do  not  have  the  values  of  derivatives  /'(())./'(/,).  However,  they 
can  he  approximated  by  liniU  differences  using  only  the  values  of  f(.r).  To  preserve  the 
correct  order  ot  accuracy  for  a  cubic  spline  approximation,  we  suggest  using  the  following 
approximations 

m  =  }X>/(A-/>)  +  <7(/'5)  (2..d7) 

"  r  =o 

/U)  =  -j  +  o(hs). 

1  k= 0 

where  h  >  0  and  p  >  d.  For  p  =  d,  if  we  take 

1  1 

co  =  .  f’i  =  d 

(> 
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Ci  ~  3 


then  .s  =  3  in  (2.37),  and  thus.  equ'>,'~>n  (2.33)  is  satisfied  within  an  error  of  0(ft3).  Corre¬ 
spondingly.  the  coefficients  o,  '  ^  <  1  for  L ,.jf(.r)  heroine 


f>i  =  v'kf(kh),  o2  =  f(0) 
k=0 
V 

O;)  =  -  fi/<  L  ~  M  ).  =  /(  L  ) 


(2.38) 


where 


1  ■  i 

?7T,0'-o)' 


2-'+1  h 


2J  +  l  h 


\  <  k  <  }>■ 


Now  we  have  /  ( ./• )  -  L,tJ/(.r)  £  //,;(/)  and  the  decomposition  (2.30)  can  he  applied  for 
it.  Therefore,  we  can  find  an  approximation  f.(.r)  for  any  function  /(.r)  £  //"(/)  as  close  as 
possible,  provided  that  j  is  large  enough,  in  the  form  of 


/, (.r)  —  Is.;  /  i-  /o  -f  fju  +  </ 1  •  ■  •  -f  !'h 


(2.30) 


where  /0(.r)  £  Vo.  </,  6  H).  0  <  /  <  j. 


3  Discrete  Wavelet  Transform  (DWT) 

In  this  section,  we  will  introduce  a  fast  Discrete  Wavelet  Transform  ( D\\  1)  which  maps 
discrete  sample  values  of  a  function  to  its  wavelet  interpolate  expansions.  Such  expansion 
with  the  wavelet  decomposition  will  enable  us  to  compute  an  approximation  of  the  derivatives 
of  the  function. 

Interpolant  Operator  Iv  0  in  V'0 

Consider  any  function  f(.r)  £  //,j ( / )  and  denote  the  interior  knots  for  l<’>  by 

.r[-1)  =  A-.  k  =!.•■•./.-  1  (3.1) 

and  the  values  of  /(.r)  on  {x[.  by 


k  =  1 .  •  •  ■  .  L  -  1 . 


I'll. 

■  ruhir  interpolant  Iv  ,,/(./■)  ot  dat a  {  /'  11 } 

i  an  he 

1  as  follows. 

1  -  1 

I I'd. /  (  r  )  -  c- 

.  i  O’.i ./'  i  r  ^  r 

.  '.'’iii  l  -r 

)  4-  (  i  .  jf.'ii  |  /  -  ./■  ! 

Chi) 

A.’  =  t  ) 

and 

;  I vii/l-c)  interpolates  <lata  l[ 

/.  -  1. 

namelv 

Iin./  ( 

/,  - 

1 

-  1. 

cut 

Let  B  he  i  lie  translorm  matrix  between  f1 

1-11  __ 

.  ci  - 1 ) 

t./| 

.  /  j  _'| 1 1  and  the 

coetln  lent 

,  ,  T  ■ 

r  - 

1  r_  |  j )  .la'. 

j-i  - 1 1  _ 

Be 

id.  hi 

win 

are 

(  h  h 

i  i  i 

\ 

i.  :*  a 

i  2 

t, 

y 

B  = 

j  j 

t ,  ^ 

1 

V 

1 

t . 

j  | 

l  £  t 

111  order  I  o  oht  a  i  II  the  Coe|  1 1 

i  lent  s  c., .  —  1 

/.•  v' 

/.  d  in 

1  d ..'5 1 .  \>. e  ha',  i  ■  1 1 

1  M  'he  the 

!  i  iad  i  a  m>n  a  I  system  Lh">  I  which  im « > I \ i  '•/.  i  operat  L ins. 

Interpolation  Operator  I,, .J  in  11  , 

Similarly,  we  ran  define  t  he  interpohit  ion  operator  I,, ,  /  ( ./• )  in  11  , .  j  >  0  for  any  ft  met  ion 
/’(.;•)  in  / /,f  ( / ) .  I*  or  this  purpose.  we  choose  the  following  interposition  points  in  /. 

..(./)  _  k  +  1  ■’ 


■>.i 


-1  <  k 


:s<i) 


where'  Uj  ~  DitnW,  =  2; I. . 

It  ran  he  easily  cheeked  that  the  interpolation  points  { ,r\  **}  tor  U,  in  (-5.1)  and  {.rjd1} 
for  ll  >  0  in  (•{.(>)  satisfy  a  “point-wise  orthogonality'  condition. 

Point  Orthogonality  of  {.r^}  for  j  >  i.  — !  <  k  <  n}  —  2 . 

CjA-r^)  =  1 

plik(.r{l'])  =  0.  -i  <  f  <  ji,  -  2  if  /  >  0;  1  <  /  <  L  —  1  if  /  =  —  I .  (-*1.7) 
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This  orthogonality  conditm.i  will  he  crucial  in  obtaining  a  fast  Discrete  Wavelet  Trans¬ 
form  (DWT). 

Tin'  interpolation  lW)f(.r)  of  a  function  f(.r)  £  in  U  J%j  >  U  can  be  expressed  as  a 

linear  combination  of  k  =  —  1 .  •  ■  •  ,  Vj  —  2.  namely 

nl~'* 

lu,f(.r)=  E  /**>(*) 


I -o/(4J))  = /(4- ^  -1  <  A-  <  m,  -  2. 

If  we  denote  Mj  as  the  th  order  matrix  which  relates  f(j>  —  (/,,_],•••.  fj.n}--i)T  anc> 

f(j)  =  (/( J-j), '  '  •  -  /Un,U))T<  then 

f(j)  =  M,f(j)  (3.9) 


where 


Mj  = 


The  solution  of  the  coefficients  -1  <  k  <nj -2  again  involves  solving  tridiagonal  system 
(3.9)  which  costs  (5 n3)  operations. 

Now  let  us  assume  that  the  values  of  a  function  f(x)  €  H'o{I)  are  given  on  all  the  inter¬ 
polation  points  {x[j)}  defined  in  (3.1)  and  (3.6),  we  intend  to  find  the  wavelet  interpolation 
Vjf(x)  C  V'o  ®  W0  0  W\  ■  ■  ■  ©  Wj  for  J  >  0,  i.e. 

L-4 

Vjf{x)  =  0  E  f-1,k<t>k(x)  +  —  L  ~  x) 

k= 0 

+  El  53  A**m*)] 

k~-\ 


/-.(*)  +  £/>(•»■) 


(3.10) 
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where 


11  j  - 

/-I  (  J')  —  Ivo/U) 

€  1  u  -  fj  (  r 

)  =  Y  €  Wj.j  >  0, 

A  —  —  1 

=  /(4~n)- 

1  <  k  <  L  -  1 

Pj/lT’l 

=  /(4J)). 

J  >0.-1  <  A-  <  tij  -  2. 

(3.11) 

Let  us  denote  f  =  (f*  f|0h  •  •  • ,  f<J>)T  the  values  of  f(.r)  on  all  interpolation  points,  i.e. 

f'-"  =  {/(dr"»£;,‘. 
fljl  =  {/(T'lKiT  i>  o. 

and  f  =  (fl-1\ f*0),  •  •  •  .  f(J))T  the  wavelet  coefficients  in  the  expansion  (3.10) 

f'-"  = 

t,j>  =  {./Uii'i'T  j>o. 

The  following  algorithm  provides  a  recursive  way  to  compute  all  the  wavelet  coefficients 
f.  and  also  the  wavelet  expansion  (3.10)  can  he  expanded  as  needed  to  include  higher  level 
wavelet  spaces  il  ).J  +  1  <  j  <  .1'  l>y  adding  only  terms  from  tin'  higher  wavelel  spaces,  i.e. 

H ■  -  H>. 

DWT  transform 
f  —  f 

This  direction  of  transform  is  straightforward  !»y  evaluating  the  expansion  (3.10)  al  all 
the  collocation  points  {.rjd'j.j  >  —1  to  obtain  f.  The  "Point-wise  Orthogonality"  (0.7) 
of  the  interpolation  points  and  the  compactness  ol  .sup/n.y  t(.c)  can  lit'  list'd  to  reduce  the 
number  of  evaluations. 

Sumbir  of  Opi  rations 

Let  .V  he  the  total  number  of  collocat  ion  points  anti  A  =  (  /,  —  1  )  +  £',_„  » ,  =  2',+t  /.  —  1 .  In 
the  evaluation  of  Vjf(.r^),  values  of  t.’(.r)  ami  o(.r)  at  tlyatlic  points  jj.d  <  k  <  2'  L.  j  >  (I 
are  needed  and  they  can  be  computed  once  for  all  for  future  use. 


Recalling  (3.10)  and  the  '‘orthogonality  condition”  (3.7)  of  the  interpolation  points,  we 
have 

7V(4f")  =  /- >(4~").  i<i<£-‘ 

which  needs  7(L  —  1)  (flops). 

For  each  0  <  j  <  J\  to  compute  Vjf{x^]),  -1  <  k  <  n3  -2,  it  needs  (5j  +  12)*Tij  (flops). 
Thus,  it  takes  7(L  -  1)  +  £/=0(5j  +  W)ni  =  2J+l  L{bJ  +  7)  +  5  L  -  7  <  7NlogN  {flops)  to 
compute  the  vector  f. 

f  — 4f 

Recalling  that  f  =  (f(-1), f(0),  •  •  • ,  f(J))T,  we  proceed  to  construction  of  Vjf{x)  in  the 
following  steps. 

Step  1 
Define 

1-4 

/_ i(x)  =  IvofM)  =  +  X]  f-\,k<t>k{x)  +  -  x), 

k= 0 

so  interpolates  /(x)  at  the  interpolation  points  x[  1\— 1  <  k  <  L  -  1,  namely 

=  M'");  (3-12) 

Step  2 
Define 

no  -2 

AM  =  I,(f|0>  -  (I m/)10')  =  E  fojtoA-r)  (3.13) 

(  =  -l 

where  (Ivo/)(0)  =  {^vof(^k))Yk0=-\ 

As  a  result  of  the  ‘point-wise  orthogonality”  conditions  (3.7)  of  the  interpolation  points, 
we  have  tko,i{x\  ^ )  =  0,  —  1  <  /  <  no  —  2,  \  <  k  <  L  —  l.  thus 

fo(x[-'])  =  0,  \  <  k  <  L  —  l. 
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So  we  have 


/-.<4"")  +  /o(4~")  =  /-i(4'")  =  W(4'")  =  /(4"") 

/-i(40,)  +  /o(40))  =  W(4“’)  +  (C  -  (W)fi  =  ff  = /(401).  <-J-<4) 

Equation  (3.14)  implies  that  function  /_i(x)  +  /o(x)  actually  interpolates  /(x)  on  both 
interpolation  points  {x* for  Vo  and  the  interpolation  points  {x^0*}  for  Wo- 

Step  3 

Generally,  we  define  for  1  <  j  <  J 

!M  =  I„,(f w  -  CP,-,/)1’1)  (3.15) 

n]  ~  2 

=  E  hk'l’iAx)-  (3.16) 

fc=-1 

where  (Py-i/)**  =  Pj-i/(**  >),  -I  <  k  <  iij  -  2. 

Again,  as  in  step  2  we  can  verify  that  function  /-i(x)  +  /0(x)  +  •  •  •  +  /y(x)  interpolates 
function  /(x)  on  all  interpolation  points  {x*-1*},  •  •  •  {x^}.  Especially,  for  j  —  J  we  have 

Vjf(x)  ~  /_ i(x)  +  /o(x)H - h  fj{x),  which  will  satisfy  the  required  interpolation  condition 

(3.11). 

Number  of  Operations. 

For  j  —  —1,  the  number  of  operations  to  invert  (3.9)  using  Thomas  algorithm  to  obtain 
is  5 L{flops).  For  0  <  j  <  J,  the  cost  of  computing  the  coefficients  in  /_,( x )  = 
M/(>)  -  (^-i/)(j))  =  E-i  <k<nj-2  fj,k^j,k(x)  consists  of  three  parts:  (1)  evaluation  of 
('Py-i/)^*  =  {Vj-if(x^)}  —  (5j  +  7)iij(flops)-,  (2)  calculating  the  difference  /b)  —  (Vj-if)^ 
-  rij(flops);  (3)  inverting  the  matrix  My  in  (3.9)  -  5nj(flop.s),  totaling  (5 j  4-  13)«y( flops). 
So  the  total  cost  of  finding  f  —  5L  +  Ey=o (6j  +  13)«y  =  E/=o(^i  +  13)2J  L  <  6N  log  N  where 
again  N  =  2J+l  L  —  1 . 

Now  let  us  go  back  to  (3.10)  to  see  the  meaning  of  the  wavelet  coefficient  {f3,k}  in  the 
finite  wavelet  decomposition  of  space  Hq(I)  for  function  /(. r).  For  this  purpose,  we  first  take5 
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a  look  at  the  wavelet  coefficient  in  the  finite  wavelet  decomposition  of  space  !*(/).  i.e.  in 

the  decomposition  \)  =  If,  H  o  -rU’i - r  h  j-  The  orthogonality  her*'  is  in  the  sense  of  Ll 

norm  not  of  H^{I)  norm.  For  simplicity  .  we  .-till  use  the  notation  o(j~)  and  i. ■  ( .r )  to  denote 
the  scaling  function  and  the  wavelet  function  for  this  decomposition,  while  keeping  in  mind 
that  they  have  different  definitions  from  those  of  //q(/)  ■  Then,  we  can  write  the  wavelet 
coefficient  in  the  finite  decomposition  as, 

fjk  =  Ji  f(x)v'jk(T)dx 

where  {r/>h.}  is  the  dual  wavelet  basis  of  {v,*.}  in  Wr  i.e.{^*fe}  is  such  a  basis  of  that 


From  Theorem  A,  we  ran  claim  that  the  absolute  value  of  the  wavelet  coefficient  \fjk\ 
depends  upon  the  local  regularity  of  f(x)  in  the  neighborhood  of  the  abscissa  2 ~Jk\  More 
precisely,  if  2~Jk  G  (a,  6),  the  decay  of  I/,*]  depends  upon  the  Lipschitz  regularity  of  f(x) 
over  the  interval  [a,  6],  as  the  resolution  2J  increases.  This  property  of  the  wavelet  coeffi¬ 
cients  allow  us  to  detect  the  location  of  the  singularity  of  the  function  and,  then,  provide  a 
general  knowledge  of  the  distribution  of  wavelet  basis  functions  whose  coefficients  are  larger 
in  magnitude  than  a  given  threshold.  The  detail  can  be  referred  to  [11].  In  the  framework 
of  L2,  the  wavelet  function  4>(x)  has  at  least  1  vanishing  moment.  Hence  the  property  of  the 
wavelet  coefficients  mentioned  above  is  always  valid. 

Now  we  return  to  the  wavelet  coefficients  {fjk}  in  (3.10).  It  can  be  easily  checked  that 

7*  W 

ip(u)  =  -(2  -  cosu;)(— — -)4e  2  . 

’  4 

Hence  ip(0)  =  which  implies  that  i/>  has  no  vanishing  moment  at  all  (see  Figure  3)  .  Since 
the  wavelet  decomposition  we  considered  here  is  in  the  space  Hl(I),  therefore,  the  decay 
property  for  the  wavelet  coefficients  fjk  ought  to  be  related  to  the  vanishing  moments  of  the 
second  derivative  of  rp(x)  (seee  Figure  4),  uot  to  those  of  i/»(x).  We  shall  illustrate  this  more 
precisely. 

Let  {?/>*fc}  be  the  dual  basis  of  {4>jk}  in  Wj  (recalling  that  the  space  we  consider  here  is 
Hi)  .  It  can  be  proven  that  for  ty**.,  (3.17)  and  (3.18)  still  hold.  Then  we  have 

/ifc  =  lf"(*Wk)"(x)dx.  (3.20) 

Notice  that  spline  wavelets  i(>(x)  and  1 l’b{x)  defined  by  (2.19)  and  (2.20)  are  continuous 
and  their  second  derivatives  have  2  vanishing  moments.  Then  applying  Theorem  A  to  the 
second  derivatives  of  ?/>(•£ )  and  ?/>(,(. r)  (namely,  by  taking  <l(x)  in  Theorem  A  to  be  j/’"(.r)  and 
?/’"( j~)  respectively),  we  can  prove  the  following. 

Lemma  4  Let  0  <  o  <  1  and  /  G  //y(  l).  If  the  second  derivative  of  the  function  /  is  Holder 
continuous  with  exponent  o,  at  .r0  G  /.  i.e. 

1/  (x)  —  f  (.r0)|  <  (  |.r  —  ./'«)[ ",  x  G  ( -c 0  —  b.  ,r()  +  b)  C  / 
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for  some  8  >  0,  then  for  any  k  6  Z,j  €  Z+  such  that  2  Jk  £  (x0  —  8/ 2,  Jo  +  A/2). 


|/,.t|  =  0(2-‘Q+1)j)!  as  j  — »  oo. 


Proof.  We  have  ,  since  if  and  t/’6  defined  by  (2.19)  and  (2.20)  are  continuous  and  their 
second  derivatives  have  2  vanishing  moments,  by  Theorem  A. 

l/>*l  =  I  <  /-  ifk  >  I  <  5Z  l°wl  I  <  />  tfj/  >  I  +  Y  l«*i)  I  I  <  /'  Vji  > 

2~}  l£(xo-6,xo+6)  2~U£(xo-6,i0+t>) 

<  Y  A'A|A-,|(9(2-(ft+,)j)  +  £  A'Ai'_a,C  =  <9(2_(o,+1)j) 

'2-J/€(r0-6,i0+«)  |/-*r|>2j| 

where  C  in  the  first  but  last  equation  is  a  constant  which  depends  on  the  second  derivative 
of  f{x). 


Lemma  4  implies  that  the  wavelet  coefficients  j  >  0,  still  reflect  the  singularity  of  the 
function  to  be  approximated.  In  practice,  when  we  solve  PDE’s  using  collocation  methods, 
we  often  use  the  values  of  the  functions,  not  their  derivatives.  Therefore,  in  order  to  use  the 
wavelet  coefficients  to  adjust  the  choice  of  wavelet  basis  functions,  we  have  to  establish  a 
relation  between  the  magnitude  of  the  wavelet  coefficients  >  0  and  f{x).  Let  us  first 

state  the  following  result  on  the  inverse  of  tridiagonal  matrix  from  [15]. 

Lemma  5  Let  A  be  a  nxn  tridiagonal  matrix  with  elements  «2,  «3,  •  ■  •  ,an  on  the  subdiagonal, 
b\,  bz,  •  •  • ,  bn  on  the  diagonal  and  q,  c3,  •  •  • ,  cn  on  the  superdiagonal,  where  c,,  c,  ^  0.  Define 
the  two  sequence  {um},  {t>m}  as  follows: 


Uo  —  1>  ktm  ( U tii _ I  Hjji—2  T  An— i  i  )  in  P  2 

c,n 


(3.22) 


^u+i  ffi  1,  vm  ( An + 1 Um ]  T  oni_|_2 1 )  in  T  ii  1  (3.23) 

^m  + 1 

where  ax  and  cn+1  are  arbitrary  nonzero  constants.  Then  A~x  —  (cv,-j)  is  given  by 


«.-j  = 


_ J-T  £Jl 

a  i  vo  ^  *  njc 

_  p 

fi  i  vo  11  a* 


(3.24) 


Corollary.  Let  A/_,  be  the  interpolation  matrix  in  (3.9),  then  we  have  rhe  lollowing 


estimates  on  Mj  1  =  («,._,), 


K,l  < 


K 


a1- 


where  K  =  1.1726  and  a  =  7  4-  \/l92  =  13.928. 
We  delay  the  proof  of  (3.25)  to  the  Appendix. 


Theorem  3  Let  f(x)  €  H$(0,L)  and  M  =  7naxi\f(x)\  and  IW]j{x)  be  its  interpolation  in 
Wj  defined  in  (3.8)  and  if  for  e  >  0,  —  1  <  k\  <  k2  <  n3  —  2 


\f{x[j))\  <  t  for  kx  <k<  k2. 

then  define 

Lj(*)=  £  MiA*),  (3-26) 

where  /  =  /(e)  =  min(^,  We  have 

|i«,,/(x)  -  I„,/(x)|  <  C(M)<  (3.27) 


where  C(M )  =  +  M),  K  =  1.1726  and  a  =  7  +  \/T92  =  13.928. 

Proof.  From  (3.9),  we  have 

f(j)  = 


where  f(j>  =  (/,  ,  •  ■  • ,  /,,„,_2)T,  fU)  =  (/(*- }),  •  •  • , f(*nJ-a))T,  thus 

/>.*  =  £>«/(*«>.  -1  -2. 

1=1 

So  we  have 

(3.28) 

For  any  given  e  >  0,  we  take  ^  =  min(rij/ 2,  —  log  e/ log  a).  For  k  €  [A,'i  +  £ ,  —  f],  using 

(3.28)  we  have 


<  /'[  E 


l/(x‘,_’!)i+  E 


!*-«!<' 


cdA:  — 


alk  — 
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<  A'<  £ 


-n - -  +  MK  T  — 

\k-t\<(alk~l\  \k%talk 


<  2A'f[l +  (-)  +  ■•• +  (-/]  + 2.1/ A'[(-)/+1  + 

ft  ft  a 

-  2Ac  — .  -  +  2.1/  A'(  -  )'+1 - - 


l-(^) 


a 


1  -  (i) 

v  O  ' 


<  2  A'  f  — +  2Af  K  e  -  1 


ft  —  1 


ft  —  1 


=  <7'f 


+  (-)’ 
ft 


where  C'  —  ^-(ft  +  A/). 
Finally,  we  have 


\Lj(x)  -  I„,/(X)|  =  |  £  h^,Ax) I 

<  H  l/j.fc|fe(*)l  <  C'c  l^j, *.•(*)) 

fc€[*-i  +^,*:2  — e]  k€[ki+(,k2-(] 

Note  that  in  the  las*  summation,  only  three  terms  will  be  nonezero  for  any  fixed  x,  so 
we  have 

\L}f(x)-lW}f(x)\<^C'(  =  Ce 

where  C  =  ~\(ot  -f  M).  This  concludes  the  proof  of  the  Theorem. 


□ 


Remark.  As  a  consequence  of  Theorem  3,  the  coefficients  fhk  of  the  wavelet  interpolation 
operator  IWjf(x)  can  be  ignored  if  x[^  E  [xj^+/,  xj^_f]  where  the  function  f(x)  is  less  than 
some  given  error  tolerance  e.  This  procedure  will  only  result  in  an  error  of  0(e).  For 
e  =  1 0— 10 ,  T  =  9,  e  =  1 0— 8,  T  =  7.  In  the  wavelet  interpolation  expansion  (3.10),  1^  is 
used  to  interpolate  the  difference  between  a  lower  level  interpolation  Vj--\  f{x)  and  /(x), 
i.e.  Vj-\f(x)  —  f{x).  Thus,  the  situation  mentioned  here  will  occur  in  larger  region  of  the 
solution  domain  as  j  becomes  larger,  avoiding  adding  unnecessary  expansion  terms  t/>j>(j). 
This  fact  will  be  used  in  the  later  section  to  achieve  adaptivity  for  the  solution  of  PDE’s. 
The  idea  of  decomposing  numerical  approximations  into  different  scales  has  been  previously 
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used  successfully  in  the  shook  wave  computations  with  uniform  high  order  spe<  iral  methods, 
where  ENO  finite  difference  methods  anti  spectral  methods  art'  combined  to  resolve  the 
shocks  ami  the  high  frequencey  components  in  the  solution,  respectively  [17]. 

We  conclude  this  section  with  the  following  result  which  shows  how  to  use  wavelet  coef¬ 
ficient  to  estimate  the  data  interpolated  by  Iu,  . 

Theorem  4  Let  IWjf(x)  and  f(x)  as  in  Theorem  4.  Anti  if  for  <  >  0,  —  1  <  k\  <  k-2  <  ii}  —  2 

|/i.fc|  <  f  for  £i  <  k  <  k2, 

then 

|/(4J))i  <  3f  for  *i  +  3  <  k  <  k2  -  3.  (3.29) 

Proof.  The  proof  follows  from  the  definition  of  IWjf(x). 

□ 


4  Derivative  Matrix  V 


The  operation  of  differentiation  of  functions  ,  which  are  given  in  terms  of  the  wavelet  ex¬ 
pansion  of  (2.39),  can  be  represented  by  a  finite  dimension  matrix  V.  Such  matrix  has  been 
investigated  in  [16]  for  wavelet  approximation  based  on  Daubechie’s  compactly  supported 
wavelets  for  periodic  functions.  The  properties  of  matrix  V,  especially  of  its  eigenvalues, 
affect  very  much  the  efficiency  and  stability  of  the  numerical  methods  for  the  solution  of 
PDE’s  to  be  discussed  in  the  next  section. 

We  consider  the  derivative  matrix  which  approximates  the  first  differential  operator 

Cu  =  uT  (4.1) 


with  the  boundary  condition 


u(L)  =  0. 


(4.2) 
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Because  of  the  inultiresolution  structure  of  spaces  lj,i.e.  Ij  C  f;+i  and  \  o  -r  IV  u  -r  ■  •  •  T 
W’j  =  Vj  y i .  We  can  rewrite  the  wavelet  interpolation  vj(x)  of  (2.39)  for  function  u(x)  as  a 
linear  combination  of  Ibj+if(x)  and  basis  in  \  J+ 1 '  namely 

V-4 

Uj(x)  =  IbtJ+iu{x)  +  ii_i(pbj+1{x)  +  ^2  uk<j>j+i,k(x)  +  uL'-sOb,j+\(L  -  x)  (4.3) 

b—Q 

where  L'  =  2'/+1  L  and  lbj+\u(x)  is  defined  in  (2.33). 

With  the  transformation  £  =  2'/+1x,  equation  (4.3)  becomes 

L'-4 

uj{£)  :  =  Uj(x)  =  hMt)  +  “-1  <Pb{£)  +  UkMO  +  UL'-3Qb(L'  -  £)•  I4--1) 

*.=0 

Using  the  notations 


u'  =  (w(l), «(2),  ■  ■  ■  ,u(L'  —  1))T  6  RL-\ 


u  =  (u(0), (u,)T,u(Z/))T  €  Ri  +1 

U  =  (u_;  ,  UOt  '  '  -  t  Ul>~3)r  6 

and  equation  (3.5),  we  have 

u  =  B-,(u'-u„)  (4.5) 


where  vector  Uf,  is  defined  by 


ufe  =  (Itiou(l),0,---,0,I6.ou(L'-  1))T  €  R^-1 


and 


Ifc,ou(l)  =  ^(co,c'1,c'2,C3,0,---,0)u/ =  7iu',  7i  €  RL  1 
0 

lbfiu(L’  -  1)  =  ~(0,  •  •  •  ,0,  -C3,  -c'2,  -c\,-Cq)  =  7-2u\  72  €  R 


L'-\ 


Therefore, 


u  =  B-’CI- 


7i 

0 


0 

72 


)u'  =  B-Tu' 


(4.6) 
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where  I  is  the  ( L'  —  1 )  x  (L1  —  1 )  identity  matrix. 

To  obtain  approximation  to  the  derivatives  of  n(£),  we  differentiate  equation  (1.1)  with 
respect  to  £  and  evalute  at  ^  =  k\  0  <  k  <  L' .  i.e. 

L'-a 

=  (I*.ow)'(6)  +  )  -  ul--aO[(L'  -  6) 

A:=0 

=  w'|(^Jt)  +  u'ldk)  o  <  k  <  /.'  (1.7) 

where  u\(^)  denotes  the  first  term  in  t he  first  equation  am  I  «';(&)  the  rest .  Recalling  the 
definition  of  Iy0u(£)  it>  (2.T5)  and  coefficients  o*.  in  (2.TS)  with  //  =  1  and  j  =  .1  +  \ ,  l.  =  /.', 
we  have 
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H  B  Tu' 


(0.  //B“T.0)u 


Finally,  combining  equations  (4.8)  and  (-1.9),  we  have 


Dili  A  0) 
DiU  A  1) 

Di<i A  i') 


=  P'u 


(4.10) 


where  the  derivative  matrix  V  is  defined  bv 


V  =  A  +  (0.  H B-T.O). 


Converting  to  tlie  j’-derivatives.  we  have 


DrUj(x I ) 

DTUAXL') 


_  •)  1 


/>€*O(0) 


=  2/+,r>'u 


(4-12) 


\  DruAXL')  /  \  DiUj(L')  ) 

where  .r,  =  jhrr-  0  <  i  <  L'. 

Let  V  be  the  upper  left  U  x  V  submatrix  of  V ,  jt+tV  will  be  the  wavelet  derivative 
matrix  to  differential  operator  (4.1)  with  boundary  condition  (4.2).  namely.  V  will  maps  the 
function  values  u(*o)>  ?<(-rih  •  •  •  ,  u(.r /,._,)  to  its  derivatives  u'(xo).  »'(* ,).  •  •  •  . 


u'(.r0) 

“'(•Cl) 


-  2j+]T> 


u(.r0) 

*4*1 ) 


(4.13) 


\  u'(*Z.'-l)  /  \  u(*Z.'-l)  / 

In  Figure  15,  we  plot  the  eigenvalues  of  T>  for  L  =  8,  7  =  0, 1,2,3  which  corresponds  to 
W  =  8,  16,32,64.  The  eigenvalues  come  in  conjugate  pairs  with  two  pure  real  eigenvalues. 
The  real  part  of  all  the  eigenvalues  are  negative  and  except  one  eigenvalues,  all  the  rest  are 
close  the  imaginary  axis. 

5  Adaptive  Wavelet  Collocation  Methods  for  PDE’s 

In  this  section  we  consider  a  collocation  method  based  on  the  DVV  1  transform  given  in 
Section  3  for  time  dependent  PDE’s.  Let  u  =  u{r,t)  be  the  solution  of  the  following  initial 


value  problem 


Hi  +  fr (  »  )  =  Mj-j-  +  </(  u  ).  -r  t  [o,  L].t  >  0 

*dlU)  =  </„(')  r 

»(/../)  =</i(/) 

1  a(j',U)  =/(./■) 

where  only  Direchlet  boundary  eonditioiis  are  considered,  however,  the  methods  presented 
here  ran  also  he  modified  to  treat  Von  Neuman  type  or  Kobin  type  boundary  conditions. 

We  use  the  idea  of  method  of  lines  where  only  the  spatial  derivative  is  discretized  by  the 
wavelet  decomposition.  The  numerical  solution  iij(.r.t)  will  be  represented  by  an  unique 

decomposition  in  Vu  :  Hu': . r  Wj.J  >  0,  namely 

L-  i 

Uj(.i\t)  -  I hJu(x,t)+  «_,  _,(/)0fc(.r)  +  +  u_i,/,-:dO'.W,( L  -  x) 

k=i) 

J  n,-2 

j  = 0 k=- 1 

J 

=  u.i(x)  +  ^2uj(x)  (5.2) 

j=0 

where  I i,,ju(x,t)  given  in  (2.3d)  consists  the  nonhomogenuity  of  u(xj)  on  both  boundaries, 
and  the  coefficients  uJtk(t)  are  all  functions  of  t.  Using  the  DWT  transform,  we  can  also 
identify  the  numerical  solution  uj(xJ)  by  its  point  values  on  all  collocation  (previously 
named  interpolation  )  points,  i.e.  { ,r [J ' }  in  (3.1)  and  (3.6),  we  put  all  these  values  in  vector 
u  =  u(/.),  i.e. 

U  =  u  ( / )  =  (u(-u.u(,,).---.u(J))T 

where  u*7'  =  {»(.r[J*./)}.  \  <  k  <  L  —  \  for  j  =  —  1:  —  1  <  k  <  n,  —  2.  for  j  >  0. 

To  solve  for  the  unknown  solution  vector  u ( / ) .  we  collocate  the  PDK  (5.1)  on  all  colloca¬ 
tion  points,  then  we  have  the  following  semi-discretized  wavele;  collocation  method. 

Semi-Discretized  Wavelet  Collocation  Methods 


»./,  +  fA'lj)  =  «./„  +  </(«/)|ra.H..  I  <  j  £  I 

n ,/(()./)  =«/(,(/) 

ii.i(L.I)  =  //,(/) 

».i(r  =  =  /(•''  = 

where  1  <  k  <  L  —  1  for  /  =  —  1;  —  1  <  k  <  n t  —  2,  for  j  >  0 
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Equation  (5.3)  involves  total  (2J+1  —  1  )L  +  2  unknowns  in  u  two  of  which  will  be  de¬ 
termined  by  the  boundary  conditions  and  the  rest  are  the  solutions  of  the  ODE  system 
subject  to  their  initial  conditions.  In  order  to  implement  the  time  marching  scheme  for  the 
ODE's  system  (for  example  Runge-Kutta  type  time  integrator),  we  have  to  compute  the 
derivative  term  in  (5.3)  fx{uj(x[^))  anti  ujrx(x[J'>)  in  an  efficient  way.  Let  us  only  discuss 
the  first  derivative  which  involves  the  computation  of  the  nonlinear  function  f{uj(x,  t.)).  For 
this  purpose  we  first  find  a  similar  wavelet  decomposition  as  (5.2)  for  f{uj).  For  a  general 
nonlinear  function  /(u),  this  can  be  done  quite  straightforward  using  the  DWT  transform 
in  Section  3. 

Computation  of  fx{{x[} *)  =  fx(uj(x[n) 

Step  1  Given  u  =  (u(_1),  u(0),  •  ■  • ,  u(-7))T,  compute  f(j)  =  {/(ulJ))},  j  >  -1  and  define 

f  =  (f(-1)if(°))...if(J))T. 


Step  2  Compute  the  wavelet  interpolation  expansion  using  DWT  transform  for  f, 

L- 4 

fj(x,t)  —  Ifc.j/  +  +  f-l,L-l{t)<t>b{L  —  x) 

k=0 

j= 0  *.•=-! 

Step  3  Differentiate  (5.4)  and  evaluate  at  all  collocation  points  {x  > -i, 


(5.4) 


L- 4 


AJ)\ 


k= 0 


+  fiAt)ViAxk])\- 

1=0  /=- 1 

Cost  of  Computing  the  Derivatives. 


For  each  single  collocation  point,  it  takes  7  +  5 (.7  +  1)  =  5.7  +  12( flops)  to  compute 
Therefore,  the  total  cost  of  computing  all  derivatives  is  (5.7  +  1 2 ) A  <  5N  log  N. 


Again,  and  o'(x)  at  the  dydir  points  0  <  k  <  2 1  L  ran  he  precomputed  once  and 

for  ail. 

Assuming  that  Euler  forward  is  used  to  discretize  the  time  derivative  in  (5.3),  we  obtain 
a  fully  discretized  wavelet  collocation  method. 

Fully  discretized  Wavelet  Collocation  Method 


=  Ujn  +  A<[-/,(u3)  +  UnJxx  +  </(Uj)]|T=rO),  -1  <  J  <  j 

«3(0)  =g0(n 
unj(L)  =p,(t“) 

“°(*  =  x{kj))  =  f{x  =  x[j)) 


(5.5) 


where  1  <  k  <  L  —  1  for  j  =  —  1;  —  1  <  k  <  rij  —  2,  for  j  >  0  and  tn  =  nAt  is  the  time 
station  and  At  is  the  time  step. 


Adaptive  Choice  of  Collocation  Points 

In  equations  (5.2)  and  (5.4),  uj(x)  and  J(uj(x))  are  expressed  using  the  full  set  of 
collocation  points  {x^}.  As  discussed  in  the  remark  after  Theorem  3  of  Section  3,  most 
of  the  wavelet  expansion  coefficients  uhk  for  large  j  can  be  ignored  within  a  given  tolerance 
f.  So  we  can  dynamically  adjust  the  number  and  locations  of  the  collocation  points  used  in 
the  wavelet  expansions,  thus  reducing  significantly  the  cost  of  the  scheme  while  providing 
enough  resolution  in  the  regions  where  solution  varies  much.  We  can  achieve  this  adaptivity 
in  the  following  two  ways. 

Deleting  Collocation  Points 

Let  e  >  0  be  a  prescribed  tolerance  and  j  >  0,  i  —  f(e)  =  min(^,  —  log ef  logo). 

Step  1.  First  we  locate  the  range  for  the  index  k, 

(fci»*i),”-,(*m,C),m  =  m(j,  f)  (5.6) 

such  that 

|uj,fc|  <  f,  k\  <  k  <  /',  i  =  1,- •  •  ,m.  (5.7) 
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Step  2.  Following  Theorem  3  and  4,  we  can  ignore  Uj^  in  Uj(x)  in  (5.2)  for  k,  <  k  < 
l,,i  —  I ,  •  •  • ,  m,  k ,  =  k'  +  (  +  3,  lt  —  l[  —  (  —  3,  namely  we  redefine  Uj( x)  as 


Uj(x)  :=  Uj.k'I’jA*) 

-1  <k<’h 

where  Kf  =  Ui 

Step  3.  The  new  collocation  points  and  unknowns  will  be 

=  1,  -,T  -  1  if;  =  -1;A;  €  {-!,-•  ■  ,n}  -2}\^. 


Increasing  Level  of  Wavelet  Space. 

Let  f  >  0  again  be  some  prescribed  tolerance,  and  if 

>  f 


(5.8) 


where  subscript  n  indicating  the  solution  at  time  t  =  f\  then  we  can  increase  the  number 
of  wavelet  spaces  Ws  in  the  expansion  for  the  numerical  solution  uj(x)  in  (5.2),  say,  up  to 

Step  1  At  t  =  tn  if  condition  (5.8)  is  satisfied,  let  J'  >  J  and  define  a  new  solution  vector 

u’j,  :=  (u<-,),u(°>,-,uW,u,;+,),-,u<J'))T 

where  for  J  +  1  <  ;  <  J',  u ^  =  {'Pju(x^)}^”2, . 

Step  2  Use  Uj,  on  the  right  hand  side  of  scheme  (5.5)  to  advance  the  solution  to  time 
step  <n+1  and  obtain  solution  Uj*1.  Then,  UjA(x)  =  Vj'Unjtx  €  Vo  0  W0  0  •  •  •  0  Wj>  will  be 
the  new  numerical  solution  which  yields  better  approximation  to  the  exact  solution  of  (5.1). 

6  Numerical  Results 


CPU  Performance  of  DWT  transform 

The  theoretical  estimates  of  operations  for  performing  the  DWT  transform  in  both  di¬ 
rection  and  the  computation  of  derivatives  at  all  collocation  points  are  0(N  log  N)  where  N 
is  the  total  number  of  terms  in  the  wavelet  expansion  (3.10). 
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We  take  the  function  in  (6.1)  and  define  its  wavelet  interpolation  expansion  (3.10)  for 
L  =  10,  J  =  2, 3,  •••,9.  t  he  total  number  of  terms  (or  collocation  points)  .V  =  2J+1  L  —  1 
are  between  79  and  10240.  In  Figure  5  we  plot  the  CPI  time  for  the  performance  of  DWT 
back  and  forth  in  both  directions  (‘o’  in  the  Figure)  and  the  computations  of  derivatives  on 
all  collocation  points  ('  +  ’  in  the  Figure).  Also  drawn  in  the  Figure  is  a  straight  line  which 
indicates  a  almost  linear  growth  of  the  CPU  timing  up  to  10 k  points. 

Adaptive  Approximation  of  Wavelet  Interpolation  Expansion 


We  consider  function - 


/(*)  = 


h\{x  +  1,0.3) 
0 

h\(x  +  0.5, 6) 
0 


if  -  1  <  x  <  -0.7 
if  —  0.7  <  x  <  —0.5  —  6 
if  -  0.5  -6  <  x  <  -0.5  +  6 
if  —  0.5  +  6  <  x  <  0 


sin(57rx)/ii(x  —  0.25,0,25)  if  0  <  x  <  0.5 
MtP)  if  0.5  <  x  <  1 


(6.1) 


where  6  =  0.01  and  hx(x,a)  is  an  exponential  hat  function  and  h2(x)  is  a  step-like  function 
and  they  are  defined  as 

_  /  exp(-^~ •)  if  |x|  <  a 


h\(x) 


0 


otherwise 


(6.2) 


and 


hz(x) 


0 


dt 

1 


if  x  <  0 
if  0  <  x  <  1 
otherwise 


(6.3) 


First  we  construct  the  full  wavelet  interpolation  expansion  (3.10)  Vjf{x)  for  J  =  6,  L  = 


40,  the  total  number  of  wavelet  functions  (or  the  collocation  points  A’’  )  Ar  +  4  =  (2 J+x  L  — 
1)  +  4  =  2J+1  L  +  3  =  5123  (including  four  boundary  functions  in  Ibjf(x)).  In  Figure  6,  on 


the  top  we  plot  the  /(x)  (solid  line)  and  Vjf(x)  at  non-interpolation  points,  at  the  bottom 


we  have  the  absolute  error  in  logarithm  scale.  In  Figure  7,  we  plot  the  components  /0  €  Uo 


and  <7j(x)  G  Wj, 0  <  j  <  6  in  Vjf{x)  =  li,jf[x)  4-  /o  +  go  +  •  •  •  +  (Jj.  We  can  see  that  only 
higher  frequency  part  is  retained  in  higher  wavelet  spaces  Wj  (notice  that  the  scales  varies 
in  different  pictures). 
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Then,  we  use  the  procedure  at  the  end  of  Section  5  to  filter  out  those  coefficients,  thus 
deleting  the  corresponding  collocation  points,  f]k  which  are  less  than  t  in  magnitude.  In 
Figure  8,  we  take  e  =  10-5  and  the  number  of  wavelet  functions  fJtk  reduced  to  289  with  the 
accuracy  of  the  approximation  (bottom  curves)  within  order  of  t.  In  Figure  9,  we  plot  the 
solution  at  the  remaining  interpolation  points  and  the  expected  clustering  of  the  interpolation 
points  is  seen  at  locations  where  the  function  changes  more  dramatically.  In  Figure  10,  we 
plot  the  magnitude  of  the  wavelet  coefficients  fj,k,j  >  —1  one  level  above  another.  High 
density  of  the  wavelet  coefficients  reflects  the  existence  of  high  gradients  of  the  approximated 
function.  In  Figure  11,  we  take  e  =  10-4  and  the  number  of  wavelet  functions  /y*  reduced 
to  206  with  the  accuracy  of  the  approximation  (bottom  curves)  within  order  of  e. 

Linear  Hyperbolic  PDE's 

We  consider  the  IVB  problem  of  linear  hyperbolic  partial  differential  equation 

ut  +  ux  =  0,  0  <  x  <  1 

•  u(0,t)  =  C  (6.4) 

.  u(*,0)  =  Mm) 

where  b  =  0.05  and  h2(x)  is  defined  in  (6.3). 

We  apply  the  collocation  method  with  adaptive  choice  of  the  collocation  points  L  = 
20,  J  =  4.  Second  order  Runge-Kutta  method  is  usd  for  the  time  derivative.  With  every 
10  iterations  we  change  the  number  and  locations  of  the  collocation  points  according  to  the 
criteria  proposed  at  the  end  of  Section  5.  The  cut-off  tolerance  t  =  I0-5.  The  number  of 
collocation  points  involved  fluctuates  ar  >imd  200  in  contrast  to  the  full  set  collocation  count 
which  is  610  in  this  case.  In  Figure  12.  we  plot  the  numerical  solution  ('+")  against  the  exact 
solution  ('o’)  at  time  /  =  0.1.  In  Figure  13,  we  plot  the  errors  in  logarithm  scale  (notice  the 
y-seale  starts  at  -2  which  corresponds  to  an  error  of  10~2).  Again,  we  see  the  automatically 
clustering  of  the  collocation  points. 

In  viscid  Burger  Equation 
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Finally,  we  consider  the  IVB  problem  of  nonlinear  hyperbolic  partial  differential  equation 

'«*  +  (£),  =  0,  - 1  <  -r  <  2 

<  u(0, t)  =  given  (6.5) 

.  u(x.O)  =  f(x) 

where 

,,  \  _  f  7.sin(7rx)  if  —  1  <  .r  <  1 
|  0  otherwise 

In  this  case,  we  take  Z,  =  10,  J  =  6.  Second  order  Runge-Kutta  method  is  used  for 
the  time  derivative.  With  every  10  iterations  we  change  the  number  and  locations  of  the 
collocation  points  according  to  the  criteria  proposed  at  the  end  of  Section  5.  The  number 
of  collocation  points  involved  fluctuates  around  100  in  contrast  to  the  full  set  collocation 
count  which  is  1280  in  this  case.  The  cut-off  tolerance  c  =  10~5.  In  Figure  14,  we  plot  the 
numerical  solutions  at  time  t  =  0.05, r  1.  The  numerical  scheme  automatically  puts  more 
collocation  points  near  the  high  gradient  (x=0)  and  the  derivative  discontinuity  (x=l). 

7  Conclusion 

In  this  paper,  we  have  constructed  a  fast  Discrete  Wavelet  Transform  (DWT)  which  enables 
us  to  study  collocation  methods  for  nonlinear  PDE’s.  The  adaptivity  of  wavelet  approxima¬ 
tion  is  conveniently  implemented  through  the  examination  of  the  wavelet  coefficients.  The 
preliminary  tests  on  the  solution  of  PDE’s  indicates  such  an  approach  will  be  important 
in  large  scale  computation  where  the  solution  develops  extremely  high  gradients  in  isolated 
regions,  and  uniform  mesh  is  not  practical.  Such  investigations  are  actually  being  done  for 
reacting  flows,  the  results  will  be  reported  in  a  separate  paper. 
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Appendix 


Proof  of  (3.25).  The  proof  is  a  straightforward  application  of  Lemma  5.  For  A/,  in 
(3.9),  we  have  (a2,  a3,  •  •  • ,  «„)  =  -75, ' ' '  *  ~n)<  (6i >  ■  ■  •  -  K)  =  (1, 1,  - • ,  1)  and 

(c-2,c^,  •  ■  •  ,c„)  =  (  — Tj,  —  75,  ■  •  •  ,  -75,-75)  where  n  =  ?ij  —  2J  L.  Therefore,  the  sequence 
{u„.}  in  (3.22)  satisfies  the  following  relations, 

2534 

u0  =  0,  u,  =  1,  Ui  =  14,  u3  =  — (A.l) 

and  for  4  <  m  <  n 

um  =  cm  (  um_2  +  14um_,)  (A. 2) 

where  cm  =  if  rn  =  n,  cm  =  1  otherwise. 

Recursive  relation  (A. 2)  is  a  finite  difference  of  order  2  whose  general  solution  is  of  the 
following  form 

«,n=c,„(ci«m-3  +  c2/?m-3)  (A. 3) 

where  a  =  7  +  \/l  92/2,  ()  =  7  —  x/ 192/2  are  the  two  distinct  roots  of  the  quadratic  equation 

x2  —  14a:  +1=0, 

and  constant  cj  and  c2  are  chosen  so  equation  (A. 3)  is  valid  for  in  —  2,3. 

Therefore, 


wo  =  0,  =  1 

=  cm(/i,o”‘-3  +  /c2/3"*-3),  2  <m<n  (A. 4) 

where  /i,  =  ^(^  “  W  >  0,/ia  =  ^(14or  ~  tt)  >  0. 

Similarly,  we  can  show  that 

0,  vn  —  1,  (A.  5) 

=  c?1_m+1(/i,o”_2_m  +M2/^n"2_m),  for  I  <  77i  <  n  -  1  (A. 6) 
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and 

t’0  =  --(Mn‘4  +  Mn~4). 

a  1 

where  <5,  =  >  0, 62  =  HMslhi  <  0. 

Finally,  following  (3.24)  in  Lemma  5,  we  have  the  following  estimates  on  the 


My 


Denote  e;,  1  <  j  <  n  as 


f  1 


13 


I  11 

V  13 


if  J  =  1 
if  j  =  2 

if  3  <  j  <  n  —  1 
if  j  =  n. 


Case  1:  i  <  j  and  1  <  j  <  n  —  1,  i  =  1 


(^icvn  2  J  2  J) 

a,J  “  _e*c"-J+1  (6ia«-4  +  Mn"4) 

So  we  have 

I  <  14an-2-J(/*i  +M2*n-2-J) 

l^ol  -  13  an-4(^_^n-4) 

<  14  o(^i  -f  ^2/ z)  1 

~  13  (<$1  -  <52)  aJ_1 


where  2  =  ^  <  1  and  K%  =  1.1666. 

Of  — 


Case  2:  t  <  j  and  1  <  j  <  n  —  1 , 2  <  1  <  n 

_  (micv*-3  +  /x2/3-3)(/i1on-2-J  +  H20n-2~*) 
Qi'>~  e,Cn-J+lC‘  (^a”-4  +  <52/?"-4) 


Case  3:  i  <  j  and  j  =  n,  i  =  1 


_ 1 _ 

e,(6ian-4  +  fa,8n-4) 


Case  4:  i  <  j  and  j  =  n,  2  <  i  <  n 


-  3  +  3) 


inverse  of 


(A.7) 


(A- 8) 


(A.9) 


(A. 10) 
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Case  5:  i  >  j  and  j  =  1,  1  <  i  <  «  —  1 


f  }Cn~x  +  \ 


/i,o 

(6,on_4  +  62(ln~4) 


Case  6 :  i  >  j  and  j  =  1 ,  *  =  n 


_ 1 _ 

^(^o"-4  +  8i(in~4 ) 


Case  7:  i  >  j  and  2  <  j  <  n  —  1 , 1  <  i  <  n  —  \ 

(/HOJ-3  + 


®i,j  —  ejcjcn-i+1 


(8^an~4  +  62/1n-4) 


Case  8:  i  >  j  and  2  <  j  <  n  —  1 ,  i  =  n 

_  +  toft3"3) 

a,J~  e}{6,a"~4  +  Mn_4)‘ 

For  Cases  2  -  8,  we  can  similarly  obtain 


K;l  < 


Kj 


2  <  i  <  8 


where  A'2  =  1.1726,  A'3  =  1.1607,  A%  =  1.1666,  A's  =  1.1666,  K6  =  1.1607,  K7  = 
A's  =  1.1666. 


Finally,  if  we  choose  K  =  1.1726,  then 


oh-’l  ’ 


1  <  *,7  <  n- 


This  concludes  the  proof. 


(A. 11) 


(A. 12) 


(A. 13) 


(A. 14) 


1.1722  and 


(A- 15) 


□ 
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Figure  1. Interior  scaling  functions  <f>(x )  (top)  and  boundary  scaling  function  &(x). 
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Figure  2.Interior  wavelet  functions  0(i)  (top)  and  boundary  wavelet  function  Tpi ,(x). 
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Figure  3,  Fourier  Transformations  of  ip(x) 
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Figure  4.  Fourier  Transformations  of  i>"{x) 
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Figure  5.  CPU  timing  for  Performing  DWT  (both  directions)  transformation  (‘o’)  and 
computation  of  derivatives  ('+’),  solid  line  -  Linear  fitting. 
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Figure  6.  Wavelet  approximation  of  function  (6.1)  with  L  =  40,  J  = 
solution  (solid  line)  and  approximation  (‘o’);  Bottom  -  absolute  error  in 
Total  number  of  is  5123. 


0.8 


6.  Top  -  Exact 
logarithm  scale. 
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Figure  9.  Close  up  of  top  part  of  Figure  8,  numerical  solutions  (‘+’)  at  remaining  collo¬ 
cation  points  against  exact  solutions  (‘o’). 
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Figure  11.  Same  as  Figure  8,  but  with  e  =  10-5.  Total  number  of  fjjt  left  is  206. 
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Figure  12.  Adaptive  collocation  solution  of  linear  PDE  (6.4)  at  f  =  0.1  with  L  =  20,  J  —  4 
and  error  tolerance  c  =  10~4.  Total  number  of  collocation  points  is  around  200. 
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Figure  13.  Error  of  numerical  solution  in  Figure  12  in  logarithm  scale. 
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Figure  14.  Adaptive  collocation  solution  of  non-linear  PDE  (6.5)  at  t  =  0.05,0.1  with 
L  =  10,  J  =  6  and  error  tolerance  e  =  10-5.  Total  number  of  collocation  points  is  around 
100. 
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Imaginary  (*2*\J)  Imaginary  (*2*\J) 


Real  (*2**J)  Real  (*2**J) 

Figure  15  Eigenvalues  for  the  first  derivative  V  for  L  =  8,  J  =  0, 1, 2, 3  whose  sizes  are 
2 J L  —  8, 16, 32,  and  64,  respectively. 
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